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ABSTRACT 

The stellar initial mass function (IMF), along with the star formation rate, is one of the fundamen- 
tal properties that any theory of star formation must explain. An interesting feature of the IMF is 
that it appears to be remarkably universal across a wide range of environments. Particularly, there 
appears to be little variation in either the characteristic mass of the IMF or its high-mass tail between 
clusters with different metallicities. Previous attempts to understand this apparent independence of 
metallicity have not accounted for radiation feedback from high-mass protostars, which can dominate 
the energy balance of the gas in star-forming regions. We extend this work, showing that the frag- 
mentation of molecular gas should depend only weakly on the amount of dust present, even when 
the primary heating source is radiation from massive protostars. First, we report a series of core 
collapse simulations using the ORION AMR code that systematically vary the dust opacity and show 
explicitly that this has little effect on the temperature or fragmentation of the gas. Then, we provide 
an analytic argument for why the IMF varies so little in observed star clusters, even as the metallicity 
varies by a factor of 100. 

Subject headings: ISM: clouds — radiative transfer — stars: formation — stars: luminosity function, 
mass function — turbulence 



1. INTRODUCTION 

The stellar initial mass function - the distribution of 
stellar masses at birth - appears to be remarkably con- 
stant accross a wide range of star- forming environments. 
While there is s ome evidence that the IMF is different in 
giant elliptical (van Dokku m fc Conro y 2010; Treu et al. 
20101 and dwarf galaxies ( Lee et al.||2009|), m ou r own 
galaxy the IMF is practically uni versal (Kroupa [2002 



galaxy the imi 1 is practical ly uni versal (ivroupa zvvz 
Chabricr 2003; Bastian et al.||2010 ). In particular , while 



zero-metailicity stars likely had a m uch different mass 
distribution than p resent-day stars (Abel et al. 2002 
|Bromm et al.| [2002), the IMFs observed tor population . 
and 11 stars appear to v ary little with the m etallicity of 
the star-forming region. jBastian et al. (2010) review the 
literature on this topic and find no evidence for a sys- 
temic dependence of the IMF on metallicity. For exam- 
ple, in the Milky Way disk, extreme outer galaxy clusters 
with metallicities of ~ 0.1 solar have the same IMF as 



nearby clusters with approximately solar metallicity (Ya- 
sui et al.||2008|). Likewise, globular clusters, which have 



metallicities that are often 10-100 times lower than the 
solar value, have the same IMF as regions of present-day 
star fo rmation once dynamical evolut ion is taken into ac- 



al 



2000 2010). 
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galaxy, |Schmalzl et al.| (|2008p and|Sirianni et al 
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find no difference between the Milky Way disk IMF and 
that of the Small Magellenic Cloud, which has a metal- 
licity of ~ 0.2 solar. In short, the IMF appears to be 
insensitive to metallicity across a wide range galactic 
environments, and understanding why is an important 
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challenge for theories of star formation. 

One pote ntial explanation was offered by |Elmegreen| 



et al. ( 2008 1 . They focus on the characteristic mass of the 
IMF plateau, which should be related to the Jeans mass 
Mj oc T 3 / 2 p -1 / 2 of the region. They calculate Mj for a 
prestellar core under the assumption that T is set by the 
balance between molecular cooling and heating from gas- 
dust coupling, finding that for a given dust temperature, 
the gas temperature depends only weakly on metallicity. 
However, this skirts an important question: what if the 
dust temperature itself is a strong function of metallic- 
ity? In regions where stars are already forming, radiative 
feedback should dominate the dust's energy balance even 
before the onset of nuclear burning, especially when the 
stars are massive. This is not a minor effect. Most stars 



are born in clusters (Lada & Lada 2003) , and the clus 
ter mass function is dN/dM oc ( [Zhang fc Fall|l999 



Fall et al. 2009 Chandar et al. 20101), which implies equa 
mass per logarithmic bin. Thus, 2/3 to 3/4 of stars are 
born in clusters of 1000 Af Q or more, and essentially all 
of these clusters are expected to contain at least one early 
O star. Hence, the maj ority of stars do n ot form in iso- 
lated environments like Elmcgr een et al.| consider, but 
form rather in the presence of a massive star that will 
affect the gas temperature distribution. If this heating 
were strongly metallicity dependent, then Mj would be 
as well. 

Other work has considered the effect of protostellar 
feedb ack on the f orm and apparent universality of the 
IMF. Bate (2009) argued that the lack of variation in 
the IMF is the result of self-regulating feedback from 
radiating protostars, but did not explain why this effect 
should be independent of metallicity when the material 
aroun d the stars is optically thick. Additionally, while 
[Bate] includes radiative feedback, he underestimates the 
luminosity by a factor of 20 and the temperature in the 
core by more than a factor of 2 (Of fner et al.||200"9 ). He 
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also does not form any massive stars, which if present 
would dominate the heating rate. 

High-mass stars are known to strongly affect the tem- 
perature and fragmentation of molecular gas. Heating 
from their extremely high luminosities, which can reach 
10 5 times the solar value, is capable of raising Mj sig- 
nificantly and preventing further fragmentation. Obser- 
vations using ground-based interferometers reveal that a 
number of ~ lOOAf massive cores appear to remain sin- 
gle, compact objects when observed at 1700 AU resolu- 
tion, even though the cores initially c ontained roughly 
10 3 Jeans masses worth of mater ial (|Bontemps et ah 



20101. Similarly, Longmore et al. (201T) show observa- 
tions of a young, massive star-tormmg region within the 
IRAS 18032-2137 complex that suggest a ty pical frag- 
ment size o f > XMq. Numerical simulations (Krumholz 
et al.||2010l hereafter KCKMIO) find that this ettect is 



most pronounced in regions of high surface density, where 
the accretion rates are high and the core is effective at 
trappin g radiation in the optically th ick regions near pro- 



tostars. 



Krumholz & McKee ( 2008 1 argued for an effeo 



tive threshold surface density above which a region will 
fragment into a few massive objects rather than a few 



sponds so well to the core mass function (e.g., 


Alves et al. 


2007 


Enoch et al. 


2008 


Beuther & Schilke 


2004) even 



tor massive 10 M ) cores, which based on isothermal 
assumptions should fragment into many small objects in- 
stead of a few massive ones. 

However, the strength of this effect could in principle 
depend on the opacity of the dust, since that determines 
the matter-radiation coupling, and hence on the mctal- 
licity of the region. It is not obvious that the above ef- 
fect would work at all in globular clusters, which have far 
fewer metals, and presumably far less dust, than present- 
day Milky Way star-forming regions. The purpose of 
this paper is to gauge the importance of metallicity to 
the fragmentation of star-forming molecular gas at small 
scales, where the collapse is highly non-isothermal due to 
radiative feedback. To do this, we have conducted a se- 
ries of core collapse simulations where we vary the dust 
opacity and show that it makes little difference to the 
core's temperature or its fragmentation. We also pro- 
vide a simple analytic argum ent based on the work of 
Chakrabarti & McKee (2005) for why the temperature 
profiles of prestellar cores should depend only weakly on 
the dust opacity. This work sheds light on why the IMF 
should be so similar in regions with different metallicity, 
even when radiation feedback from massive stars on the 
gas cannot be ignored. 

The outline of the paper is as follows. In Section [2] we 
describe our simulation setup, including the initial con- 
ditions and numerical methods used. In Section [3] we 
show the results of our simulations, which demonstrate 
how fragmentation is insensitive to metallicity variations 
of a factor of 20. In Section [4] we apply the work of 



Chakrabarti & McKee ( 2005[ ) to show that the tem- 
perature profile of a dusty, centrally-heated core should 
depend only weakly on its metal content. Finally, we 
present our conclusions in Section [5] 

2. SIMULATION SETUP 
2.1. Numerical Techniques 



To perform our simulations, w e have used ORION 
(Truelove et al. 1998 Klein 1999), a parallel, adaptive 
mesh, radiation-hydrodynamics code for astrophysical 
applications. The equations and methods used are al- 
most identical to those in KCKMIO, so we describe them 
only briefly here, emphasizing the differences. In short, 
ORION tracks the four conserved quantities of the com- 
bined gas-radiation fluid: the density p, the momen- 
tum per unit volume pv, the non-gravitational gas en- 
ergy density pe, and the radiation energy density E. It 
updates these quantities conservatively, including the ef- 
fects of gravity and diffuse radiation, but not magnetic 
fields or ionizing radiation. The radiation transport is 
solved in the flux-limited diffu sion approximation u sing 
the mixed-frame formulation ( Mihalas fc Klein|1982 ), re- 
taining terms of order v/c. For a full description of the 
equations solved by the radiation module as well as accu- 
racy tests of the mixed- frame treatment, see |KrumhoIz] 
etaL] p307b| and |Shestakov fc Offner|p008| . 

In addition to the gas and radiation, the simulation do- 
main also includes star particles, which are placed on fine 
level grids whenever the Jeans density in a cell exceeds 
a critical value (see below). Star particles interact with 
the gas by accreting mass fro m the simulation volum e ac- 



cording to the algorithms in Krumholz et al. ( 2004 ) and 
emitting radiation according to the protos tellar model 



described in the appendices of Offner et al. (2009). The 
full set of equations solved for the fluid is 



d_ 

Ft 



p=-V • (pv) - ^2 MiW(x - x t ), 



dt 



(pv) = -V ■ (pvv) - VP - pVc/) - WE 
~^2PiW(x- Xi), 



(1) 



(2) 



d_ 
Ft 



(pe) = -V • [(pe + P)v] - pv ■ - K 0P p(4irB - cE) 



-A 2— -1 v VE- V£w(as-a%), (3) 



dt 



E = V 



cX 



VE ] + kopp(4ttB - cE) 



korp / 
M2 K °^-Av-VE-V-( 3 F R2 vE 



Hon 

^LiW(x- Xi), 



(4) 



where «or and kqp are the Planck and Rosseland mean 
opacities computed in the co-moving frame of the fluid 
and B = caaTg / (4tt) is the Planck function. The flux 
limiter A and R2, which is related to the Eddington fac- 
tor, are dimensionless numbers that enter into our ap- 
pr oximation of the radiati ve transfer. For details, refer 



to [Krumholz et al. (2007b) and the references therein. 

In the above equations, the sums are taken over all 
the particles in the simulation. Li is the luminosity of 
star i , while Mj, p i: and £j are the rates at which mass, 
momentum, and energy are transferred from gas to stars. 
W(x — Xi) is a kernel that distributes the transfer over 
a radius of 4 fine level cells around the star. The star 
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particles are updated according to 



dt 

d_ 
dt 

d 
dt 



-Mi = Mi. 



Xj - 



Pi 



Mi 



-MiVcj) + Pi , 



where <p is the gravitational potential given by 



= -AttG 



p + 2_ J M i 8(x- as,-) 



(5) 
(6) 
(7) 

(8) 



For the purpose of computing the gravitational force on 
a sink particle, we have used a Plummer Law with a 
smoothing length of 2 fine level cells (14 AU in physical 
units) for the gas-sink force and a smoo thing length of 1 
fine le vel cell for the sink-sink forces; see |Krumholz et ah] 
( 2004 ) for details. Note that we do not include the effects 



of protostellar outflows; for discussion of high-mass sta r 
formation with outflows, see Cunningham et al. (2011). 
We adopt a polytropic equation of state: 



P = 



^k = (7-l)p(e- 



(9) 



where T g is the gas temperature, /i = 2.33 is the mean 
molecular weight for molecular gas with cosmic abun- 
dances, and 7 is the ratio of specific heats. We have 
taken 7 = 5/3, appropriate for a gas of molecular hydro- 
gen with the rotational levels frozen out, but this choice 
is essentially irrelevant because T g is determined primar- 
ily by radiative effects. 

The above is all identical to KCKM10. The only differ- 
ences between the numerical schemes employed in those 
simulations and ours are: 

1. We have used the dust opa city model described in 
Cunningham et all ( 1201 1[ ), based on the work of 
Semenov et al.| ( |2003p . This opacity model was cre- 
ated with stellar winds in mind and includes line 
cooling effects at high temperatures. These effects 
are irrelevant at T g < 1000 K, which includes prac- 
tically all of the gas under discussion here. 

2. We have turned off mergers between sink particles 
with masses greater than 0.05 Mq. 

3. We have modified the Plank and Rosseland mean 
dust opacities by a multiplicative constant 6 to al- 
low for dust-to-gas mass ratios other than solar. 

Item (2) requires some discussion. As mentioned 
above, each sink particle is surrounded by an accretion 
zone of 4 fine level cells from which it draws gas. In 
KCKM10, if a sink particle moved within another sink's 
accretion zone, the particles would merge regardless of 
their masses. Because the finest resolution in our simu- 
lations is ~ 7 AU, this would mean that any stars that 
ever moved within a distance of 28 AU would be com- 
bined. Because this may not be realistic, we have im- 
posed a mass limit of 0.05 Mq above which sink parti- 
cles will no longer merge. This limit is chosen to roughly 
correspond to the mass at which a protostar's core tem- 
perature becomes high enough to dissociate molecular 



hydrogen and initiate second coll apse (Masunaga et al 



1998 Masunaga & Inutsuka 2000 ) . Before this, the sinks 
are more like extended balls of gas with radii of a few 
AU than stars, so it is more likely that they will merge. 
After second collapse, the sinks represent objects that 
are much smaller, and mergers should be less likely. Al- 
though our initial conditions are also slightly different 
(see below), this choice accounts for most of the differ- 
ences between the simulations reported here and those 
in KCKM10. 

2.2. Refinement criteria 

The computational domain is a cube of side L^ ox that 
is discretized into a coarse grid of Nq cells, so that the 
resolution on the coarse grid is xo — Lt, OK /No. The AMR 
functionality of ORION automatically identifies regions 
that need more resolution and covers them with a finer 
grid. With L levels of refinement and a refinement ratio 
of 2, the resolution of the finest level is Axl is xq/2 l . 
In this work, we have chosen these parameters such that 
Axl is always ~ 7 AU. 

In generating the grids, a cell is tagged for refinement 

if 

1. It is within a distance 16Axl of the nearest star 
particle. 

2. It has a density greater than the Jeans density, 
given by 

T 2 vol 



p.i = r 



GAx\ ' 



(10) 



where c s is the soun d speed and we use J = 1/8 
( jTruelove et aL][r997l ). 



3. the local gradient in the radiation energy is greater 
than a critical value, given by 



0.15 



E 
Ax~ L 



(11) 



This procedure is repeated recursively until the final level 
is reached. Taken together, these conditions ensure that 
the regions near the star particles are always tracked with 
the highest available numerical resolution. 

2.3. Initial conditions 

Our initial conditions also follow the approach laid out 
in KCKM10, and have been chosen to resemble the struc- 
tures from which massive stars are believed to form. 
Observations of the internal structure of infrared dark 
clouds (IRDCs) using sub-mm interferometry reveal the 
presence of pe aks in the loc al density distribution termed 
massive cores (Swift 2009). They are measured to have 
masses in the range of ~ IOOMq, radii of about ^0.1 pc, 
and temperatures of about ~ 20 K. Massive cores are ob- 
served to be centrally concentrated, and unlike low-mass 
prcstcllar cores are highly turbulent, with virial ratios of 
order unity. 

To model these objects, we initialize the simulation 
volume to contain a sphere of gas with radius R, total 
mass M, surface density £ = M/irR 2 , and constant ini- 
tial temperature T C oro = 20 K. Following the theoretical 
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TABLE 1 

Simulation Parameters 



Name 


6 


M (M ) 


E (g cm- 2 ) 


R (pc) 


rr v (km s 1 ) 


iff (kyr) 


ibox (pc) 


L 


N 


Axo (AU) 


N L 


Ax L (AU) 


Solar 


1.0 


300 


2.0 


0.1 


3.59 


30.2 


0.40 


6 


192 


430 


12,288 


6.71 


0.2 Solar 


0.2 


300 


2.0 


0.1 


3.59 


30.2 


0.40 


6 


192 


430 


12,288 


6.71 


0.05 Solar 


0.05 


300 


2.0 


0.1 


3.59 


30.2 


0.40 


6 


192 


430 


12,288 


6.71 


High E 


0.2 


300 


10.0 


0.045 


5.37 


9.03 


0.18 


5 


168 


220 


5,376 


6.85 



Note. — Col. 7: linear size of computational domain. Col. 8: maximum refinement level. Col. 9: number of cells per linear dimension 
on the coarsest level. Col. 10: linear cell size on the coarsest level. Col. 11-12: same as col. 9-10, but for the finest level. 



models of McKee & Tan ( 2003 ), we give the core a power 
law density profile ol 



p(r) oc 



(12) 



where we have taken k p to be 1.5. This is consistent 
with observations of clumps of molecular gas - structures 
a few pc in size with thousands of solar masses of mate- 
rial - which reveal power-law de nsity profiles wi t h scal - 



ing exponents b e tween 1 and 2 (Bcuther et al. (2007), 



Myers (1995), Mueller et al. 1 2002 [JTandwth 



[Caselli 

higher resolution observations of individual cores, which 
find po wer law density profile s with slopes betwe en 1.5 
and 2 ( |Longmore et ;il.||201 1| [Zhang et aL"1|2009[ ). We 
stress, however, tfiat tfie artificiality of the initial condi- 
tions, in which the density initially lacks structure and 
the core is considered in isolation, is a potential source 
of uncertainty. 

The cores are placed in a cubic volume with Lbox = 
4 x i?, surrounded by a background medium with density 
Pbg = 0.01 x Pcdgoi where p e d ge is the density at the edge 
of the initial core: 



Pcdgc 



3-fc p \ M 



47T 



i? 3 ' 



(13) 



To maintain pressure balance, the temperature of the 
background gas is set to T hg = 100 x T corc = 2000 K. 
The opacity of the ambient gas is set to a tiny value to 
ensure that it does not interfere with the temperature 
structure of the collapsing core. 

To mimic the effects of turbulence, we also give the 
cores an initial Gaussian random velocity field with a 
power spectrum P(k) cx fc~ 2 . We scale each compo- 
nent of the velocity to have a one-dimensional root-mcan- 
squared dispersion equal to t he velocity at the sur face of 
a singular polytropic sphere ( McKee fc Tan|[2003 ) : 



GM 



R(2k p - 1) 



(14) 



The corresponding virial parameter a v ; r = 5<r 2 GM/ R 
( |Bertoldi fc McKee|[l992| ) is 5 for k p = 3/2, so that the 
turbulent kinetic energy is initially larger than the grav- 
itational potential energy. However, the turbulence de- 
cays as the core collapses, so the virial parameter drops. 

Starting from these conditions, we have performed a 
series of four core collapse simulations, summarized in 
Table fl] For three of the runs, we chose a surface den- 
sity ofE = 2 g cm' 2 , which is comparable to the surface 
densities observed in massive star-forming regions in the 
Galaxy. We vary the parameter 5, which represents the 
dust-to-gas mass ratio relative to that expected at so- 



lar composition. The values of 5 = 1.0, 0.2, and 0.05 
are representative of nearby star clusters, extreme outer 
galaxy star clusters, and low-metal globular clusters, re- 
spectively. We have also conducted a run at E = 10 g 
cm -2 , 5 = 0.2, which is characteristic of extragalactic su- 
per star clusters. Motivated by the above observations, 
we have chosen to set the mass M of all our cores to 
300M Q . This is higher than in KCKM10, but it will al- 
low us to form massive stars more quickly and thereby 
gauge the effect of radiative feedback at less computa- 
tional expense. Because the mass of the cores is the 
same in the E = 2 g cm" 2 and S = 10 g cm -2 runs, the 
High E run has a slightly smaller radius. In the rest of 
the paper, we will refer to these runs by the names given 
in Table [T] 

We run all the simulations out to 0.5tg, where tg is the 
gravitational free-fall time evaluated at the mean density 
p = 3M/4irR 3 : 



tg 



3vr 
32Gp' 



(15) 



By this time, the basic similarity of all the runs has been 
established. We emphasize that in general, the numerical 
methods and resolution have been held constant across 
all the runs, so that any differences between them should 
be attributable to variation in <5 or E. 

3. RESULTS 
3.1. Temperature and density structure 

Figure [T] shows the result of evolving the above cores 
out to half a global free-fall time. The leftmost panels 
show the column density in units of g cm -2 , zoomed 
out to show the entire simulation volume. The large 
scale morphology of the collapse is practically identical 
between the E = 2 g cm -2 runs, as that is set primarily 
by the magnitude of the initial velocity perturbations. 
The High E run shows a slight tendency towards more 
filamentary structure than the others, since the initial 
Mach number must be larger to pump the higher surface 
density cloud into virial equilibrium. Nonetheless, the 
differences between the High E run and the others are 
only minor at this scale. 

The middle panels again show the column density, with 
the same units and color scale as before, but zoomed in 
to show the central 5000 AU around the most massive 
object. At this scale, some differences between the runs 
are apparent. In particular, the shape of the accretion 
flow around the stars appears to be different in the three 
E = 2 g cm -2 runs. This is because radiation pressure 
can be important near the massive star(s), and its mag- 
nitude varies with the opacity. 
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le-01 1 10 100 1000 le-01 1 10 100 1000 10 100 1000 

Column Density (g / cm"2) Column Density (g / cm"2) Temperature (K) 



Fig. 1. — Projections through the simulation volumes at t = 0.5%. The left panels show the column density of the entire core, defined 
as f pdx. The middle column is also the column density, zoomed in to show the middle 5000 AU. The right column shows the column 
density-weighted temperature, f pT gaB dx/ f pdx, at the same scale. The rows, from top to bottom, show runs "Solar," "0.2 Solar," "0.05 
Solar," and "High S". Stars are represented by circles drawn on the plots, with the size of the circle corresponding to the size of the star. 
Stars with masses between 0.05Mq and IMq are the smallest, intermediate mass stars with IMq < m < 5Mq are larger, and stars with 
masses greater than 5Mq are the largest. 
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Fig. 2. — Star particle statistics as a function of time for all four 
runs. The top panel shows the total mass in stars (the top set of 
lines) and the mass of the most massive star (bottom, dotted set). 
The middle plot shows / max , the fraction of total stellar mass that 
is in the most massive star. The bottom plot is the total stellar 
luminosity. 

The rightmost panels show the column density- 
weighted gas temperature at the same 5000 AU scale. 
The most striking feature is that the temperature of the 
gas surrounding the central condensation of stars is quite 
similar in the three £ = 2 g cm~ 2 runs. Changing the 
opacity by a factor of twenty appears to have little impact 
on the temperature of the bulk of the gas. On the other 
hand, the High E run has noticeably more gas heated 
to a hundred degrees or higher. The "wall" of hot gas 
visible on the left-hand side of some of the runs is the 
hot, diffuse medium that was initially outside the core; 
it is more visible in the High £ run because the core 
is smaller, so that the 5000 AU frame captures a larger 
fraction of the volume. The stars are near the edge of 
the core at t = 0.5tg in all the runs because the random 
velocity perturbations we used happened to advect them 
that way. With a different realization of the turbulence, 
this would not necessarily happen. 

3.2. Fragmentation and Star Formation 

At the end of 0.5 free-fall times, about 40 M of gas 
has been turned into stars, or about 13% of the total 
core. At that point, all the runs have a massive star 
of ~ 10M©. In the three £ = 2.0 g cm' 2 runs, this 
star forms a binary system with a massive companion 
of ~ 6Mq. In the High £ run, the secondary is only 



3Mq , with the missing mass spread out among the other 
objects. 

A major differe nce between these simulatio ns and those 
of KCKM10 and [Cunningham et al.| ( |2011[ ), is that we 
form many more objects during the collapse. This dif- 
ference is largely due to the choice of merger criterion - 
by turning off mergers beyond a threshold mass, objects 
that would have formed a single star here form several. 
At the opening stage of the collapse, the initial velocity 
perturbations create a network of dense filaments, visible 
in the leftmost panels of Figure [T] Because the core is 
centrally concentrated, gas falls into the center, forming 
a star with about 3M Q after approximately 0.15 free-fall 
times. At this time, additional sinks begin to form in the 
gas that is falling on to the star. Some of these sinks are 
small enough that they will merge onto the central star, 
but many of them accrete enough mass that they surpass 
the 0.05M Q threshold and become "stars" in their own 
right. These objects fall into the center of the core's grav- 
ity well and begin to undergo complex N-body interac- 
tions with each other. They stay in the center for several 
orbits before being thrown out. At that time, they have 
will have grown to approximately 0.1 to 1 M Q , and will 
take away with them mass that would have merged with 
the central star(s) in the simulations without merger sup- 
pression. This results in the most massive star gr owing 
much less rapid l y after ~ 0.2tg than in KCKM10 or|Cun 



ningham et al. (2011). Thus, with mergers suppressed 
we form a massive star or binary plus a system of a few 
dozen low-mass stars, as opposed to having most of the 
stellar mass in one system. 

The fragmentation observed in these simulations is 
qualitatively different from previous simulations with 
ORION that did not suppress mergers for sinks larger 
than O.O5M . While we emphasize that the handling of 
stellar mergers is artificial in both cases, the large num- 
ber of massive st ars observed to have close companions 
( Sana et al.||2008 ) shows that merging all stars with sep- 
arations less than 28 AU is clearly unrealistic. For that 
reason, the model used here is probably more accurate 
than merging stars regardless of their masses. However, 
it probably errs in the other direction. For example, be- 
cause we necessarily soften the gravitational interactions 
of stars with gas on scales smaller than the grid scale, 
we underestimate the dynamical friction that stars ex- 
perience as they pass through one another's disks, and 
thus underestimate the rate at which this process cause s 



stars to be captured (e.g., Pfalzner et al. 



2005, 2006) 



Even if such captures do occur, our present prescription 
does not allow the stars to merge. We emphasize that 
all particle mer ger prescriptions, whether in our code o r 



in others (e.g., |Bate et al.||1995| |Federrath et al. 



are suspect because they involve subgrid models 



2010} 
or hy 



drodynamic and gravitational interactions on scales that 
are not resolved in the simulation. The absolute num- 
ber and masses of stars formed is affected by the choice 
of sub-grid model, so these results be should interpreted 
cautiously. However, because the merger criterion is the 
same for all the runs presented in this paper, any dif- 
ferences in the mass distribution of the objects formed 
between the four runs is due to changes in £ and d, not 
the merger criterion. 

The idea that fragmentation of incoming gas limits the 
mass supply for massive stars has been discussed before 
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in the literature. Most recently, Peters et al. (20101 
coined the term "fragmentation-induced starvation to 
describe the phenomenon where gas en route to a mas- 
sive star instead collapses to form a low mass star before 
arriving, depriving the massive star of that material. Be- 
cause our simulations start from turbulent initial condi- 
tions, the fragmentation happens in the dense filaments 
that feed the central star system rather than in a disk as 
in the non-turbulent simulations of Peters et. al., but the 
underlyin g idea is the same. T he competitive accretion 
models of |Bonnell et al.| ( |1997[) and |Bonnell et ah] ( |2004 1 
show a similar phenomenon. However, it is important to 
distinguish between the dynamics in our simulation and 
those in traditional competitive accretion models. In the 
absence of radiative feedback, fragmentation always pro- 
ceeds down to the th ermal Jeans mass, ~ 0.5 M Q (e.g., 
Bonnell & Bate|2006|). Only once fragments of that mass 
form do they then grow to larger masses by Bondi-Hoyle 
accretion of gas that is not initially bound to the star. 
The dynamics in our simulation are different, in that the 
radiatively heated gas in our simulations does not always 
fragment down to such small masses, and our central 
massive stars are built by a direct collapse and not by 
Bondi-Hoyle accretion. The difference beco mes apparent 
if one compares our simulations to those of |Dobbs et al.| 
(2005), who start with initial conditions very similar to 
ours, but do not include radiative feedback. In their sim- 
ulations, after O.btg they find no stars larger than ~ 1 
Mq because all the gas has fragmented to small masses, 
while in our simulations we have 10 Mq stars after a 
similar time. 

Figure [2l shows the properties of the star particles 
formed in the simulations as a function of time. To show 
all the runs on the same plot, we have normalized the 
time by tg, but in physical units the High £ simulation 
evolves about 3 times faster than the others. The top 
panel shows the total mass in stars, which is quite sim- 
ilar across all the runs. This is to be expected, because 
the star formation rate is set by the global properties of 
the flow, which are essentially independent of radiative 
effects. The middle panel shows / max , the fraction of the 
total stellar mass that is in the most massive star. There 
is a period around 0.3iff where the Solar run levels off, 
but this appears to be a temporary phenomenon. By 
QMs, / max is very similar across all the runs. Overall, 
there appears to be little difference between the runs, 
either in the total mass converted to stars or in the frac- 
tion of that mass that ends up in the most massive star. 
The bottom panel of Figure [2] shows the total luminos- 
ity of all the stars in the simulation. The value of 10 4 
Lq we find is typical o f observed massive protostars (e.g. 
Cesaroni et al.| ( |2007| ). Here, we can see a difference 
between the High run and the others - because of the 
higher accretion rates, the total luminosity is higher in 
the High £ run, although the difference decreases with 
time as the stars become more massive and more of the 
radiant output comes in the form of nuclear luminosity. 
The nuclear luminosity is never dominant, however, and 
thus the luminosity as a function of time is roughly flat 
after about 0.2tg, even though we are forming more stars 
and the stars are growing more massive. We will make 
use of the luminosity averaged after over 0.2tg to 0.5tg 
in Section 4] below. 

Figure [3 shows the cumulative mass distribution of the 



TABLE 2 

P-VALUES FROM TWO-SIDED K-S TESTS FOR THE 
SIMULATED MASS DISTRIBUTIONS. 



Run 


Solar 0.2 Solar 


0.05 Solar 


High S 


Solar 


0.1 


0.62 


0.02 


0.2 Solar 




0.16 


0.01 


0.05 Solar 






0.02 


l.Or— . ■ ■ ■ ■ ■ , 




Fig. 3. — Cumulative mass distributions from all four runs at 
t = tff0.5. /(> m) the fraction of the total stellar mass that is in 
stars with masses greater than m. 

stars in the four runs - that is, for each value of to, the 
y-axis shows the fraction of the total stellar mass than 
in is stars with masses greater than to. Visually, there 
appears to be little difference between the curves, par- 
ticularly for the £ = 2 g cm -2 runs. To test whether 
the initial mass function is indeed the same in the dif- 
ferent runs, we have performed two-sided Kolmogorov- 
Smirnov (K-S) tests between each pair of distributions. 
The results are shown in Table [2j At the 10% level, we 
cannot reject the null hypothesis that the three £ = 2 
g cm -2 have the same underlying initial mass function. 
The High £ distribution, on the other had, does seem to 

be statistically different from the others. 

In contrast t o the minor effect reported here, Krumholz 
et al.| (2007a) found that isothermal runs fragmented 
completely differently from radiative ones, and KCKM10 
found a major difference in fragmentation between runs 
with low and high surface density. To summarize, the 
differences in the fragmentation of all of our runs are mi- 
nor. To the extent that there are significant differences, 
they are due to changes in the surface density, rather 
than to changes in the metallicity. At least within the 
range of parameters considered here, metallicity appears 
to have little effect on either the temperature or the frag- 
mentation of molecular gas. 

4. DISCUSSION 

4.1. Analytic Model 

The above simulations suggest that metallicity plays 
little role in the fragmentation of star-forming gas. To 
understand why, consider a simple model system like the 
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initial conditions above: a core of gas and dust with ra- 
dius i? c , mass M, surface density £ = M/irR 2 , and a 
power law density profile p(r) oc r~ kp . We would like 
to understand what happens to the temperature of this 
core once stars have started to form, so we will place a 
point source of luminosity L in the center to represent 
the combined radiant output of the central collection of 
stars. We assume that the dust opacity follows a power 
law in the far IR regime: 



-Sk q 
-Sk 



, for 3 mm < A < 30 /im (16) 



where the subscript "0" refers to an arbitrary reference 
value and 5 is the dust-to-gas mass ratio relative to solar. 
The T in this equation is the dust temperature, which we 
assume is identical to the radiation temperature. We will 
adopt the dust opacity model of Wcingartner & Draine 
( 12001 ) , for which n = 0.27 cm 2 g 1 at A = 100 /xm and 
p = 2; however, we have verified that using the o pacities 



oflSemenov et al.j ( |2003[ ) (as used by ORION) or |Pollack 
et al. 1 T994| ) makes little difference. For reference, the 
corresponding To is 144 K. 

The emission from such a system is considered in 
Chakrabarti fe McKee| ( [2005] ) (hereafter CM2005). They 
find that even though cores do not have sharply delin- 
eated photospheres as stars do, the radiation they emit 
is still well-described by 



L = MnR 2 ch aT* h 



(17) 



where L is a constant of order unity and T c h and i? c h are 
a characteristic temperature and the radius from which 
radiation with frequency j/ c h = kT^/h has an optical 
depth of 1. These are given by 



Rr = 



(L/M)£( 4 +«/' 3 



AaL 



(3 - k p )6n 
A(k p - 1)1$ 



1 4//T 



(18) 



and 



T, 



ch 



{ 



L/M 
AaLY^p 2 



A{k p 



im 



(3 - k p )Sn Q 



(19) 



where a — 2(3 + A(k p — 1). Rather than using the ex- 
pression for L given in C M2005, we will adopt the mor e 
accurate expression from Chakrabarti fe McKee| ( |2008[ ) , 
which they report gives excellent agreement wi th results 
from the DUS TY code, based on the work of |Ivezic "fc 
Elitzur| ( p97| : 

L = 1.6i? c 01 . (20) 

A final result we will take from CM2005 is that the tem- 
perature profile in the vicinity of the photosphere is also 
well-described by a power law: 

-k T 



T(r) = T t 



ch 



\Roh) 



(21) 



We can solve the above equations simultaneously to get 
that the temperature as a function of density (or, equiv- 
alently, radius) in the core is 



T{p) = T c (p/pcdgc) 
T(r) = T c (r/R c )- kT 



k T /k p 



'L/M' 


k p -l + fik T 


"(3- 


k p )SK 


. AaL _ 




A(k p 


- 1)T$ 



4fc T -2 



(4+/3)fc T + fc p -3 



x E" 

or, specializing to k p = 3/2 and /3 = 2, 



(22) 



T 



'L/M' 


fc T + l/4 
3 


3(5 k q 


. 4:<jL _ 




. 4T o . 



2k T -l 

3~ 



We can see from this expression that the scaling of T c 
with 5 actually goes to zero for k? = 0.5. This is sur- 
prising at first, since both R c and T c h scale with S, to 
the -2/3 and -1/3 powers, respectively, for our assumed 
parameter values. However, we can understand this in 
the following way: take a point r that is outside the 
effective photosphere of the core. Since radiation with 
frequency v c h is thin here, this point's temperature will 
be determined roughly by radiative equilibrium with a 
luminosity source with effective radius R c h and temper- 
ature T c h- If we then lower 5, i? c h will decrease, since 
the core will be less effective at trapping radiation and 
photons with frequency v c h will be able to travel further 
through the core's envelope before being absorbed. This 
will make the effective emitting area smaller, which will 
tend to lower the flux at point r. However, because it 
is closer to the central heating source, the temperature 
at the new value of i? c h will be higher as well, and that 
will tend to increase the flux at r. If the temperature 
always scales as the -0.5 power of radius, then these two 
effects will cancel out exactly; the luminosity seen at r, 
L w 47reri? 2 h T c 4 h , will be the same, and the temperature 
will be the same as well. 

Is close to 0.5? CM2005 give a fitting function for 
the temperature scaling exponent kx as a function of R c : 



0.48/fcO- 005 

T - 0.02fci-° 9 
Rr 



O.lkf 5 

tir. 



(24) 



This expression shows that as long as R c > 1, meaning 
that the photospheric radius is less than half the core 
radius, kx is indeed quite close to 0.5 and depends only 
weakly on R c . This expression assumes that T c h < 300 K, 
so that most of the flux is emitted at wavelengths longer 
than 30 /im and the opacity is well approximated by a 
power law in frequency. We have verified that this con- 
dition holds under the circumstances considered in this 
paper. Since dust sublimates at ~ 1000 K, this condition 
also implies that R c h is larger than the dust destruction 
radius. 

Note that for a given S and L/M, R c is a function 
of S only. From our simulation results, we can calculate 
the light-to-mass ratio associated with E = 2 and 10 g 
cm~ 2 from the luminosities shown in the bottom panel 
of Figure [2j averaged over the period from t = 0.2iff to 
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TABLE 3 
Light-to-Mass Ratios 



Run 


L/M 


Solar 


80.7 


0.2 Solar 


69.8 


0.05 Solar 


65.4 


High S 


152.6 



Note. — Averaged over t = 0.2tff to 
t = 0.5i ff . 



2=2 g cm _i 
£ = 10 g crrr 




10°* 



10 2 



10 



Fig. 4. — Temperature scaling exponent kx as a function of <5 for 
the two different values of £ considered in the simulations. 

t = 0.5t ff . The results are shown in Table [3] The "M" 
here refers to the gas mass only, not the mass in stars. In 
the E = 2 g cm -2 runs, the light-to-mass ratio appears 
to depend weakly on d, but the effect is small and we will 
simply use the average over the three metallicities. 

In Figure |4j we plot kx as a function of 5 for the two 
valu es of E represented in our simulations from Equation 
(24). Note that kx remains relatively close to 0.5 for a 
wide range of values for 5 for both of the surface densities 
we have considered. It is therefore not surprising that the 
simulations show such a weak dependence on 6. 

Note that for most of the parameter space, the value 
of kx is actually a bit less than 0.5, meaning that the 
temperature actually scales inversely with 6. When the 
temperature falls off more slowly than kx = 0.5, the 
temperature at the smaller value of R c h associated with 
a lower 5 will actually be larger than what is required to 
keep the temperature outside of the photosphere at a con- 
stant value, so T will increase slightly. Once 5 increases 
much past solar, however, kx increases dramatically and 
T begins to rise with 5. 

4.2. Comparison to simulations 

We can also calculate the temperature profiles ex- 
pected for the cores in the above simulations. The pro- 
cedure is as follows. Using the light-to-mass ratios from 
Table[3]along with the known simulation valu es of E and 
S, we can compute kx using Equations (18) and pi] ). 

(|22l 



Then, we use this value along with Equation ([22| to get 
the temperature as function of radius. In Figure [5j we 



plot these relations and compare them to the profiles 
calculated from the simulations. To get the simulation 
profiles, we calculated the density-weighted mean tem- 
perature in spherical shells of radius r around the most 
massive star, and plotted the result versus r. The density 
weighting ensures that the hot, diffuse gas surrounding 
the cores does not interfere with the result inside the 
core, although its presence can be seen in the tempera- 
ture rise as r approaches the core radius R c . We have 
averaged the simulation profiles over all the snapshots 
from t = 0.2%, when the luminosity has roughly leveled 
off, to t = 0.5£ff , and the error bars show the standard er- 
ror over all the snapshots. The error bars are larger close 
to the central massive star because of the presence of dy- 
namically interacting stars within in central few hundred 
AU of the simulation volume. The important things to 
note are 

1. The simulation results confirm the power law na- 
ture of the temperature profile and agree closely 
with the predicted slopes 

2. The effect of varying the metallicity is quite small 
compared to the effect of changing the surface den- 
sity, both in the analytic calculation and in the 
simulations. 

The analytic expression is systematically higher than 
the simulations by roughly 10%. It does, however, agree 
with the jump in T from E = 2 g cm~ 2 to E = 10 
g cm~ 2 in quite well. The discrepancy is likely due to 
differences in the treatment of the radiative transfer be- 
tween our simulations and CM2005. ORION treats the 
radiation in a frequency-averaged gray approximation in 
terms of the Plank and Rosselan d mean opacities, whil e 
CM2005 used the DUSTY code fllvezic fc Elitlu7|[T997| , 
which includes frequency information about the photons. 
However, CM2005 also assumed a spherically symmet- 
ric, static core with no density perturbations, so it is not 
clear which result is more accurate. Whatever the case, 
both methods agree that the temperature profiles are not 
particularly sensitive to S. 

4.3. Predictions for star-forming regions 

Star-forming regions in the Milky Way and other galax- 
ies sample a large range of surface densities and metallic- 
ities. A question we can ask is: in what range of param- 
eter space is the gas temperature insensitive to changes 
in the metallicity? To answer this question, we will use 
the relationship betwee n the light-to-mass ratio an d the 
surface density given in Krumholz & McKee (2008): 



L/M pa 3.6M 2 - - 33 E!!- 67 r b j 6 



(25) 



where M 2 is M/100Af Q , E = S/l g cm" 2 , and T b ,i 
is the background core temperature divided by 10K. In 
what follows, we will take Mq, — 3 and Tb,i = 2. This ex- 
pression is calculated by assuming that the core converts 
its gas into stars at a rate of a few percent per free-fall 
time and is not accurate once massive stars have formed. 
However, it is still relevant to the early stages of collapse, 
and will allow us to gauge whether the core temperature 
becomes high enough to slow fragmentation. Note that 
the star formation efficiency assumed here is lower than 
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— Solar 

— 0.2 Solar 

— 0.05 Solar 
High Sigma 




Fig. 5. — Temperature profiles from the simulations (dots) and 
analytic theory (solid lines). To show both values of S on the 
same plot, we have normalized r by the size of the core R c . The 
simulations profiles are averaged over t = 0.2% to t = 0.5%. 



what we see in the simulations; because we have not in- 
cluded the effects of outflows or other turbulence-driving 
mechanisms, our star formation rates are higher than ob- 
served. 

Using this expression, we can eliminate the dependence 
on L in the equation for T c to get the edge temperature as 
a function of M, S and 8 only. We plot the expected tem- 
perature for a M = 3OOM0 core as a function of these pa- 
rameters in Figure [6j To charact eriz e the core by a single 
temperature, we use Equation (22) to evaluate T(p) at 
the mean density p = 3p e dge/(3- : ^p). Note that we have 
assumed the temperature is determined by protostellar 
feedback; gas with sufficiently low metallicity would be 



warm anywa y owing to the lack of coolants. [Omukai 
et al. (2010) considered this effect, and concluded that 
it should be dominant for metallicities lower than 0.01 
Zq, so we will limit our analysis to values of S greater 
than 0.01. Figure [6] includes a range of surface densities 
that span the conditions typical of star formation, from 
low-mass star-forming regions like Taurus and Perseus 
(0.1 g cm -2 ), to regions of active massive star forma- 
tion in the galaxy (1 g cm -2 ), to extra- galactic super 
star clusters (10 g cm -2 ). For all those environments, 
one would need to change 5 by at least two orders of 
magnitude to get a factor of 2 change in the tempera- 
ture. Hence, we expect there to be very little variation 
in the fragmentation of cores across environments with 
different metallicities over the range currently probed by 
observations. Note that Figure [6] shows the mean core 
temperature only when the dominant source of heating 
is protostellar feedback. In reality, cores in the upper 
left hand region of the parameter space would likely be 
hotter than the indicated 5 - 10 K due to heating from 
cosmic rays or the decay of turbulence, so the true range 
of temperature variation is even smaller than indicated 
in Figure [6] 

Unlike 0, X does make a significant difference in the 
core temperature. In regions like Taurus (E ~ 0.1 g 
cm - 2), we do not predict protostellar feedback to be 
able to heat the gas much above 10 K. Hence, collapse 
there will likely be isothermal, and may be more prone to 




Fig. 6. — Contours of the mean core temperature T(p) at the 
early stages of collapse as a function of S and 8 for a core with 
mass M = 300M Q . The turnover of the T = 30 K contour at high 
£ and 8 is due to fcy becoming large; see Equation (|24l and Figure 
II ^ 

fragmentation to produce low-mass stars. Regions like of 
higher surface density (S ~ 1 g cm -2 ), are not isother- 
mal, which may bias them towards forming high-mass 
stars. 

5. CONCLUSIONS 

We have performed a series of numerical experiments 
using ORION in which we follow the collapse of massive 
cores past the onset of star formation to the subsequent 
heating of the gas by radiative feedback, varying the dust 
opacity by a factor of twenty. We find that the opacity 
makes little difference to either the temperature or the 
fragmentation of the cores as they collapse. Our simula- 
tions consider surface densities of £ = 2 g cm -2 , char- 
acteristic of massive star-forming regions in the Milky 
Way, and £ = 10 g cm -2 , characteristic of extra-galactic 
super star clusters. 

We have also presented an analytic argument for why 
the IMF should be relatively independent of the metal- 
licity of the star-forming region, even when the heating 
is dominated by a central source as in high-mass star- 
forming cores. This helps to explain why there do not 
appear to be significant differences between the IMFs 
of disk stars and globular clusters, or between those of 
the Milky Way and the SMC, despite large differences 
in metallicity. We find that the metallicity should only 
weakly influence the IMF over a large range of star- 
forming environments. 
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